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I.  EXECUTIVE  SUMMARY 

1.  We  have  developed  a  new  hybrid  technique  embed¬ 
ding  a  quantum  description  of  point  defects  in  SiC>2 
inside  a  classical  model  of  the  infinite  solid.  Thus, 
we  incorporate  accurately  long-range  effects  (lattice 
relaxations  and  dielectric  effects).  Whilst  embed¬ 
ding  has  been  employed  previously  in  highly  ionic 
systems,  such  as  the  alkali  halides,  this  is  the  first 
successful  implementation  in  a  mixed  ionic-covalent 
system.  We  have  demonstrated  that  many  of  the 
geometrical  parameters,  such  as  bond  lengths  and 
bond  angles,  are  reproduced  consistently  in  both 
the  classical  and  quantum  regions.  Furthermore, 
important  quantum  mechanical  quantities,  such  as 
atomic  charges  and  electronic  densities  of  states, 
show  no  artifacts  due  to  embedding.  Finally,  be¬ 
cause  we  are  using  reliable  quantum-mechanical 
methods,  we  can  calculate  excited  states  with  rea¬ 
sonable  accuracy  and  so  can  study  optical  prop¬ 
erties  with  a  good  deal  of  rigor.  This  has  impli¬ 
cations  for  materials  in  electronics  and  for  optical 
fiber  materials.  While  we  have  chosen  one  specific 
material,  our  work  included  successful  creation  of 
automated  techniques  for  generating  force-field  pa¬ 
rameters  in  the  classical  region,  as  well  as  pseudo¬ 
potentials  used  on  frontier  pseudo-atoms,  so  that 
the  method  is  readily  applicable  to  emerging  mate¬ 
rials,  such  as  high-K  dielectrics. 

2.  We  have  applied  this  new  technique  to  radiation- 
induced  defects  in  both  crystalline  and  amorphous 
SiC>2  and  developed  new  insights  that  change  dra¬ 
matically  our  understanding  of  the  physics  of  these 
centers  and  hence  significantly  altered  some  of  the 
standard  models  for  defect  creation  in  radiation  en¬ 
vironments  such  as  those  found  in  space.  Specifi¬ 
cally, 

•  The  relaxation  around  the  oxygen  vacancy  is 
of  much  longer  range  than  can  be  captured  in 
current  super-cell  calculations.  This  has  im¬ 
plications  for  defect  properties  in  ultra-thin 
SiC>2  films  in  MOS  devices. 

•  We  have  measured  the  distributions  of  struc¬ 
tural  and  spectroscopic  properties  of  oxygen 
vacancy  centers  in  amorphous  silica,  that  can¬ 
not  be  demonstrated  in  other  calculations. 

•  Our  modeling  predicts  the  two  major  struc¬ 
tural  types  of  positively  charged  vacancies  ( E ' 
centers):  dimer  and  dangling  bond  centers. 


The  local  structure  of  both  types  of  centers 
depends  on  the  medium  range  structure  of  the 
surrounding  amorphous  network.  Contrary  to 
the  standard  model,  we  find  that  in  amor¬ 
phous  Si02  positively  charged  oxygen  vacan¬ 
cies  readily  distort  to  localize  spin  on  one  side 
of  the  defect  without  one  of  the  silicon  atoms 
rebonding  to  a  three  fold  oxygen.  This  implies 
that  E '  centers  should  be  much  more  common 
that  currently  thought.  For  each  type  of  de¬ 
fect  we  found  the  distribution  of  both  struc¬ 
tural  and  EPR  parameters  and  calculated  op¬ 
tical  absorption  spectra. 

•  We  have  fully  characterized  two  other  struc¬ 
tural  types  of  E'  centers:  the  dimer  and  the 
back-projected  configurations  and  calculated 
their  spectroscopic  parameters. 

•  We  formulated  the  structural  criteria  which 
favor  the  formation  of  different  types  of  center 
in  the  original  amorphous  structure  in  terms  of 
the  average  Si-0  distance  of  oxygen  ion  with 
its  two  neighboring  silicon  ions.  This  should 
allow  us  to  estimate  the  relative  concentra¬ 
tions  of  different  defect  centers  based  on  the 
mechanism  of  their  creation. 

Thus  we  have  demonstrated  that  the  method  has  ful¬ 
filled  many  of  its  promises.  It  does  give  greatly  improved 
accuracy  over  the  standard  periodic  calculations  that  em¬ 
ploy  density  functional  theory  for  the  representation  of 
many  important  experimentally  measured  quantities.  It 
has  already  generated  new  physical  understanding  of  im¬ 
portant  radiation-induced  defects,  indicating  its  poten¬ 
tial  importance  in  the  study  of  emerging  materials.  The 
results  of  this  work  are  presented  in  the  following  publi¬ 
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II.  INTRODUCTION 

Amorphous  silica,  a-SiC>2,  is  one  of  the  most  techno¬ 
logically  important  and  best  studied  amorphous  materi¬ 
als  with  a  long  history  of  defect  studies1,2.  However,  our 
understanding  of  the  effects  of  disorder  on  defect  prop¬ 
erties  and  the  relative  abundance  of  different  defect  con¬ 
figurations  in  amorphous  silica  is  only  just  starting  to 
develop.  Theoretical  modeling  has  been  instrumental  in 
developing  models  of  basic  defects  in  silica,  starting  from 
pioneering  papers  by  W.  B.  Fowler  et.al.3~7 ,  and  has  led 
to  significant  advances  in  the  field.  It  has  been  realized 
that  modeling  defects  in  amorphous  materials  requires 
considering  not  only  one  or  a  few  sites,  as  in  crystals, 
but  a  statistical  ensemble  of  structural  sites.  The  de¬ 
pendence  of  formation  energies,  structural  parameters, 
diffusion  barriers  and  spectroscopic  parameters  of  these 
defects  on  details  of  the  local  and  medium  range  envi¬ 
ronment  should  be  characterized  by  some  statistical  dis¬ 
tributions.  The  structural  disorder  can  also  lead  to  cre¬ 
ation  of  new  defect  types.  Models  of  these  in  amorphous 
materials  should  take  into  account  the  role  of  local  and 
medium  range  stress  and  strain  in  the  surrounding  net¬ 
work  in  defect  formation.  In  a  crystal,  the  defect  model 
which  comes  out  of  calculations  includes  a  unique  set  of 
its  main  fingerprints,  such  as  the  geometric  and  electronic 
structure,  spectroscopic  properties,  mechanism  of  motion 
(diffusion),  reactions  with  other  defects  and  other  param¬ 
eters  which  can  be  studied  experimentally  and  used  for 
ultimate  defect  identification. 

In  this  project,  we  have  developed  an  embedded  clus¬ 
ter  method  for  predicting  and  analyzing  defect  models  in 
a-Si02  •  The  purpose  of  this  report  is  to  critically  review 
the  scope  of  this  method  and  to  demonstrate  its  applica¬ 
bility  to  studying  the  structure  and  spectroscopic  prop¬ 
erties  of  defects  in  amorphous  silica.  Below  we  consider 
the  basic  approaches  used  in  defect  studies,  describe  the 
embedded  cluster  method  developed  within  this  project 
and  its  main  advantages  with  respect  to  other  methods, 
and  present  the  main  results  of  its  application  to  devel¬ 
oping  defect  models  in  amorphous  silica.  The  details  of 
the  computational  scheme  are  given  in  Appendix  A  and 
the  results  of  the  test  calculations  in  Appendix  B. 

III.  MODELING  TECHNIQUES 

We  start  from  briefly  discussing  some  techniques  avail¬ 
able  for  modeling  defects  in  insulators.  They  exist  in  the 
form  of  computer  codes  implementing  a  model  of  a  sys¬ 
tem  with  a  defect  and  a  method  for  calculating  the  defect 
structure  and  its  properties. 

A.  Periodic  and  cluster  approaches 

Two  models  are  often  used  in  computer  simulations 
of  point  defects  in  solids:  a  periodic  model,  and  a  finite 
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cluster  model.  They  differ  in  the  boundary  conditions 
imposed  on  a  system  wave-function.  Initially,  the  pe¬ 
riodic  model  was  developed  for  calculation  of  the  elec¬ 
tronic  structure  and  properties  of  ideal  crystals,  whereas 
the  cluster  model  was  introduced  from  molecular  calcula¬ 
tions.  However,  consequently  the  terms  ’’periodic  model” 
and  ’’cluster  model”  for  the  same  defect  appeared.  It 
has  long  been  realized  that  neither  of  them  can  be  used 
alone  if  an  accurate  picture  of  a  defect  is  to  be  estab¬ 
lished.  Therefore,  both  models  were  gradually  developed 
to  provide  accurate  descriptions  of  defect  structure  and 
processes. 


1.  Periodic  model 

The  periodic  model  is  based  on  a  definition  of  a  unit 
cell  (or  more  generally,  a  supercell),  which  can  be  period¬ 
ically  translated  to  build  an  infinite  solid  as  indicated  in 
Fig.  lb  and  lc.  There  are  no  limitations,  in  principle,  on 
the  size  and  shape  of  the  unit  cell  providing  that  the  host 
lattice  is  defined  unambiguously.  Note  that  the  same  ap¬ 
proach  is  often  used  to  build  models  of  amorphous  solids 
too.  The  periodicity  imposed  on  the  atomic  structure 
also  applies  to  the  electronic  structure:  the  charge  den¬ 
sity  is  the  same  in  every  supercell  and  is  matched  at  the 
boundary  between  two  neighboring  supercells.  In  the 
case  of  ideal  crystals  this  approach  is  technically  exact. 

To  calculate  the  properties  of  an  individual  point  de¬ 
fect  using  the  periodic  model,  a  supercell  is  defined 
which  comprises  the  defect  and  its  immediate  environ¬ 
ment  (Fig.  lb).  The  supercell  is  then  periodically  re¬ 
peated  as  shown  in  Fig.  lc.  The  defect  in  each  supercell 
interacts  with  an  infinite  number  of  similar  defects  (called 
images)  in  all  other  supercells.  Therefore,  the  application 
of  this  approach  to  individual  defects  is  adequate  only  if 
artificial  defect-image  interaction  is  small. 

One  of  the  most  attractive  features  of  the  periodic 
model  is  that  it  allows  one  to  consider  a  perfect  system 
and  one  with  a  defect  on  the  same  footing.  Most  of  the 
computer  codes  implementing  the  periodic  model  use  a 
basis  set  of  plane- waves  or  their  variations.  This  provides 
a  very  flexible,  potentially  complete  basis  set,  a  quality 
almost  inaccessible  in  cluster  calculations. 

An  obvious  feature  of  the  periodic  model,  which  can 
be  both  to  its  advantage  and  disadvantage,  is  that  ’’de¬ 
fects”  are  periodically  translated  and  interact  between 
themselves.  It  can  be  an  advantage  if  one  is  interested, 
for  example,  in  periodic  adsorption  of  molecules  at  the 
surface  and  would  like  to  vary  their  concentration.  It  is 
a  drawback  if  one  is  interested  in  the  properties  of  an 
individual  defect  in  different  charge  states.  In  particu¬ 
lar,  if  the  defect  is  charged  with  respect  to  the  host  lat¬ 
tice,  the  electrostatic  interaction  of  such  defects  is  diver¬ 
gent.  Some  special  techniques  developed  to  deal  with  this 
problem  are  discussed  in  refs. 8-10 .  Yet  another  downside 


of  the  periodic  model  is  that  the  defect-induced  distor¬ 
tion  of  the  host  lattice  can  only  be  taken  into  account 
within  the  supercell.  This  may  lead  to  inaccurate  results 
if  these  distortions  propagate  further  than  cell  bound¬ 
aries,  but  can  be  checked  by  increasing  the  supercell  size. 
Large  supercells  are  required  to  study  complex  defects 
and  to  account  for  extended  defect-induced  lattice  de¬ 
formation.  The  recently  developed  ’’so-called”  order-IV 
methods  allow  one  to  significantly  reduce  the  comput¬ 
ing  load  by  making  it  just  linearly  proportional  to  the 
number  of  particles,  IV,  in  the  system.  Supercells  includ¬ 
ing  several  thousand  atoms  can  be  considered  using  these 
methods11,12. 

Finally,  studying  excited  states  in  the  periodic  model 
represents  a  difficult  problem.  On  one  hand,  excited 
states  are  often  much  more  delocalized  than  the  system 
ground  state  and  hence  require  much  bigger  cells  for  re¬ 
liable  predictions.  On  the  other  hand,  computational 
techniques  for  calculating  excited  states  are  less  devel¬ 
oped  in  the  periodic  model  compred  to  those  available 
for  molecules  and  in  a  cluster  model. 

To  summarize,  the  periodic  model  is  ideal  if  one  is  in¬ 
terested  in  the  ground  state  properties  of  neutral  defects, 
which  weakly  perturb  the  surrounding  lattice. 


2.  Molecular  cluster  model 

An  alternative  approach  is  to  consider  a  point  defect 
and  its  immediate  environment  as  a  large  hypothetical 
molecule  usually  called  a  ’’cluster”  (see  Fig.  Id).  This 
model  is,  generally  speaking,  only  applicable  if  the  fol¬ 
lowing  two  conditions  are  satisfied:  i)  the  atomic  and 
electronic  structure  of  the  host  lattice  is  accurately  repro¬ 
duced  by  the  cluster,  and  ii)  the  defect-induced  perturba¬ 
tions  are  essentially  localized  within  the  cluster  bound¬ 
aries13. 

To  define  a  cluster  means  to  somehow  ’’cut  out”  a  frag¬ 
ment  from  an  infinitely  large  host  lattice  and  to  consider 
it  separately.  The  ’’cutting”  could  be  more  easily  made 
by  a  minimal  charge  density  surface.  This  is  relatively 
straightforward  in  atomic,  molecular  and  ionic  crystals 
where  structural  elements,  such  as  ions  or  molecules  can 
be  well  defined  and  their  electron  densities  separated.  As 
a  result,  an  integer  number  of  molecules  or  ions  can  be 
attributed  to  a  cluster.  In  the  case  of  polar  or  cova¬ 
lent  crystals  the  cutting  is  made  through  atoms  so  that 
the  cluster  should  be  terminated  by  ’’pseudo- atoms”  spe¬ 
cially  designed  to  terminate  broken  bonds  (see  review14 
for  numerous  examples). 

The  range  of  applications  for  the  cluster  model  is  of¬ 
ten  complementary  to  that  of  the  periodic  approach.  A 
cluster  model  is  a  natural  choice  if  it  is  essential  to  con¬ 
sider  the  properties  of  an  isolated  defect  in  a  crystal. 
Other  examples  include  properties  of  nano-powders,  low- 
coordinated  surface  sites  or  different  sites  in  an  amor¬ 
phous  structure,  i.e.  situations  where  supercells  are  dif¬ 
ficult  to  construct  or  too  expensive  to  consider.  One  of 
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FIG.  1:  Periodic  and  embedded  cluster  approaches  for  calculation  of  point  defects  in  crystals,  a)  A  crystal  with  a  point  defect, 
b)  Immediate  environment  of  the  defect  in  the  periodic  approach  is  defined  as  a  supercell  for  periodic  calculations,  c)  The 
supercell  is  periodically  translated  so  that  to  form  a  new  ideal  crystal,  d)  Immediate  environment  of  the  defect  in  the  embedded 
cluster  approach  is  included  in  a  quantum-mechanically  (QM)  treated  cluster,  e)  The  QM  cluster  is  embedded  in  a  classically 
treated  environment  (outer  circle). 


the  main  advantages  of  cluster  models  is  the  large  num¬ 
ber  of  quantum-chemical  methods  for  electronic  structure 
calculations  available  from  molecular  calculations,  espe¬ 
cially  for  excited  states. 

On  the  downside,  cluster  model  calculations  are  usu¬ 
ally  more  computationally  demanding  than  periodic 
ones.  In  spite  of  enormous  advances  in  computing  power, 
clusters  of  only  a  few  tens  of  atoms  can  be  calculated  rou¬ 
tinely.  This  is  because  the  computational  expense  grows 
as  n4  where  n  is  the  number  of  atomic  basis  set  functions 
typically  used  in  these  calculations. 

Using  small  molecular  clusters  for  treating  larger  sys¬ 
tems  brings  a  number  of  problems:  i)  Polarizability  of 
the  ’’surface”  ions,  especially  O2-,  is  larger  than  those  in 
the  ’’bulk”  of  the  cluster15.  A  large  ratio  of ’’surface”  to 
’’bulk”  atoms  in  a  small  cluster  may  result  in  an  overesti¬ 
mated  electronic  polarization  due  to  the  ’’surface”  atoms 
and,  in  principle,  in  an  incorrect  electronic  structure  al¬ 
together.  This  is  manifested  in  formation  of  artificial 
’’ghost”  states  in  the  band  gap.  ii)  The  electrostatic  po¬ 
tential  due  to  the  host  lattice,  which  is  particularly  im¬ 
portant  in  ionic  compounds,  is  not  adequately  treated  in 
small  clusters.  As  a  result  the  positions  of  energy  lev¬ 
els  can  be  wrong,  iii)  Defect-induced  lattice  relaxation 
can  be  considered  only  for  the  atoms  nearest  to  the  de¬ 
fect  with  atoms  at  the  cluster  surface  usually  kept  fixed. 
This  implies  that  the  lattice  response  to  the  presence  of  a 
defect  as  well  as  the  effect  of  that  response  on  the  defect 
itself  can  be  significantly  underestimated. 

In  ionic-covalent  systems,  such  as  SiC>2,  cluster  forma¬ 
tion  always  requires  breaking  chemical  bonds.  A  popular 


fix  in  this  case  is  to  passivate  broken  bonds  with  real  hy¬ 
drogen  atoms  or  with  ”pseudo”-atoms  with  parameters 
specially  fitted  to  fulfil  this  task.  A  cluster  of  S^CUHg 
including  two  tetrahedra  connected  via  a  common  oxy¬ 
gen  atom  with  all  other  oxygens  passivated  by  hydrogen 
atoms  is  a  frequent  choice.  Note  that  this  cluster  is  just 
an  artificial  construction  convenient  for  studying  the  lo¬ 
cal  properties  of  defects  in  SiCU- 

Thus  a  molecular  cluster  model  tends  to  completely 
disregard  the  effects  of  the  rest  of  the  atoms  outside  the 
cluster  and  applies  some  semi-empirical  fixes  to  the  atoms 
at  the  cluster  surface  to  diminish  their  adverse  effect. 
This  major  problem  is  addressed  and  partially  solved  in 
the  embedded  cluster  model. 


3.  Embedded  cluster  model 

The  molecular  cluster  model  can  be  improved  by  tak¬ 
ing  into  account  the  interaction  of  the  quantum  cluster 
with  the  rest  of  the  host  lattice,  the  perturbation  of  the 
lattice  by  the  defect,  and  the  effect  of  the  lattice  per¬ 
turbation  on  the  defect  itself.  This  is  achieved  by  con¬ 
structing  an  external  potential  in  which  the  cluster  is 
then  embedded.  Such  a  potential  is  called  an  embedding 
potential  and  the  model,  an  embedded  cluster  model.  De¬ 
velopment  of  an  accurate  embedding  potential  is  a  major 
challenge.  It  aims  to  ’’correct”  the  molecular  cluster  ap¬ 
proach  by  improving  the  description  of  the  interaction 
between  the  cluster  and  the  rest  of  the  system.  This  can 
be  done  at  three  levels  of  sophistication: 
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•  The  electrostatic  potential  due  to  the  part  of  the 
crystal  lattice  outside  the  cluster  is  included  in 
the  calculations  by:  i)  placing  point  charges  of 
appropriate  value  at  the  lattice  sites  (see,  for  ex¬ 
ample,16,17;  ii)  fitting  the  Madelung  potential  by 
several  point  charges  located  at  specially  selected 
points;  iii)  calculating  the  matrix  elements  due  to 
the  Madelung  potential  directly  and  adding  these 
matrix  elements  to  the  potential  energy  matrix18 . 

•  To  further  improve  the  description  of  the  cluster 
boundary,  an  interface  region  between  the  clus¬ 
ter  and  the  point  charges  is  introduced.  In¬ 
terface  atoms  can  be  approximated  by  effec¬ 
tive  core  pseudo-potentials19,20  at  the  cation  sites 
around  the  cluster  or,  more  accurately  by  effec¬ 
tive  potentials  calculated  self-consistently  specifi¬ 
cally  for  the  material  under  study  (see,  for  exam¬ 
ple,  refs.1'’21-23).  The  latter  approach  allows  one 
to  reduce  the  cluster  size  and  the  computational 
load  in  some  cases. 

•  The  defect-induced  lattice  distortion  is  accounted 
for  by  embedding  the  cluster  in  an  infinite  polariz¬ 
able  lattice  of  classical  ions  represented  using  a  shell 
model24  or  a  polarizable  ion  model25,26.  Unlike  the 
two  previous  cases,  these  ions  are  not  kept  fixed, 
but  are  allowed  to  adjust  their  positions  to  mini¬ 
mize  the  total  energy  of  the  system.  This  allows  one 
to  relax  explicitly  a  large  region  (several  hundred 
atoms)  of  the  host  lattice  near  the  defect.  More 
importantly,  embedding  into  a  polarizable  lattice 
can  be  done  in  such  a  way  that  it  also  accounts 
for  the  effect  of  the  lattice  response  on  the  cluster 
electronic  structure.  This  allows  one  to  achieve  self- 
consistency  in  the  calculation  of  the  defect-induced 
lattice  perturbation  and  the  effect  of  this  perturba¬ 
tion  on  the  defect  structure  and  properties.  This 
effect  is  similar  in  many  ways  to  a  polaron  digging 
its  potential  well  in  a  polarizable  media. 

Most  of  the  calculations  dealing  with  defects  in  amor¬ 
phous  silica  have  so  far  been  made  using  either  periodic 
or  molecular  cluster  models.  In  the  periodic  model,  the 
disorder  of  the  amorphous  structure  is  included  explic¬ 
itly,  but  within  a  relatively  small  periodically  translated 
supercell  (see,  for  example,  refs.2'-30).  In  the  molecular 
cluster  model,  a  cluster  can  capture  at  best  the  average 
local  structure  of  the  material.  A  problem  pertaining  to 
both  models  has  been  that  periodic  cells  and  clusters  typ¬ 
ically  include  up  to  72  atoms,  which  is  hardly  enough  to 
consider  the  medium-range  distortion  of  the  amorphous 
structure  induced  by  the  defect.  Another  important  issue 
concerns  the  dependence  of  the  results  on  the  quantum 
mechanical  method  (see,  for  example,  refs.31,32).  In  par¬ 
ticular,  as  has  been  demonstrated  in  refs.33,34,  the  elec¬ 
tronic  structure  of  polaron  and  exciton  systems  depends 
dramatically  on  whether  Hartree-Fock  or  Density  Func¬ 
tional  Theory  (DFT)  is  used.  Besides,  neither  of  these 


two  models  can  be  used  to  build  up  significant  statistics 
and  to  find  rare  precursor  sites.  The  major  limitations  of 
periodic  DFT  calculations  relate  to  predicting  narrower 
band  gaps  and  inability  to  predict  optical  absorption  and 
luminescence  energies  of  defects. 

The  embedded  cluster  model  is  designed  to  overcome 
the  limitations  of  both  models  and  to  allow  one  to  prop¬ 
erly  treat  the  defect-induced  distortion  of  amorphous  net¬ 
work,  build  up  defect  statistics  and  calculate  defect  spec¬ 
troscopic  properties  with  good  accuracy.  Some  examples 
of  the  development  of  this  approach  and  its  application 
to  defects  in  crystalline  oxides  are  discussed  in  refs.35-39. 
They  clearly  demonstrate  that  providing  the  system  is 
ionic,  it  can  be  well  treated  in  the  embedded  cluster 
model  discussed  above.  Several  embedding  schemes  have 
been  suggested  for  treating  ionic-covalent  systems,  such 
as  Si0231,40,41.  They  have  been  successful  in  predicting 
ground  state  properties  of  crystalline  quartz.  However, 
they  failed  to  provide  reliable  predictions  for  defect  op¬ 
tical  absorption  energies.  This  is  largely  due  to  the  for¬ 
mation  of  spurious  states  at  the  bottom  of  the  conduc¬ 
tion  band  caused  by  inaccurate  treatment  of  the  cluster 
border  within  the  local  embedding  potentials  approxima¬ 
tion31,35. 

To  summarize,  one  of  the  main  advantages  of  the  em¬ 
bedded  cluster  approach  is  that  it  is  inherently  local  by 
nature  and  can  be  easily  applied  to  bulk  crystals,  amor¬ 
phous  materials,  surfaces42-45  and  complex  interfaces46. 
However,  if  building  up  a  quantum  cluster  requires  break¬ 
ing  bonds  between  atoms,  the  Coulomb  embedding  po¬ 
tential  is  not  enough  and  a  more  sophisticated  treatment 
is  needed.  The  development  of  an  approach  to  treating 
such  systems  has  been  at  the  center  of  this  project. 


IV.  THE  EMBEDDED  CLUSTER  METHOD 
IMPLEMENTED  IN  THE  GUESS  CODE 

The  previous  generation  of  the  embedded  cluster 
method  implemented  in  the  GUESS  computer  code39 
has  been  thoroughly  described  in  several  recent  publi¬ 
cations31,32,43.  Therefore  we  will  only  briefly  outline  the 
main  features  of  the  method  and  will  focus  on  the  details 
relevant  to  studying  defects  in  a-Si02-  The  detailed  de¬ 
scription  of  the  method,  typical  calculation  setups  and 
test  calculations  are  presented  in  Appendix  A  and  B. 

In  broad  terms,  in  the  embedded  cluster  method  a 
system,  crystalline  or  amorphous,  with  a  single  point 
defect  is  divided  into  several  regions.  A  spherical  re¬ 
gion  I  centered  at  the  ’’site  of  interest”  includes:  i)  a 
quantum- mechanically  treated  cluster  (QM  cluster),  ii) 
an  interface  region  which  connects  the  QM  cluster  and 
the  (classical)  rest  of  the  solid,  and  iii)  a  classical  re¬ 
gion  that  surrounds  the  QM  cluster  and  includes  up  to 
several  hundred  atoms.  Region  I  is  surrounded  by  a  fi¬ 
nite  region  II,  also  treated  atomistically  and  containing 
up  to  several  thousand  atoms.  In  the  course  of  calcu¬ 
lations,  all  positions  of  atoms  within  region  I  are  fully 
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optimized  while  atoms  in  region  II  remain  fixed  in  the 
positions  corresponding  to  a  perfect  crystalline  or  amor¬ 
phous  structure.  Their  purpose  is  to  provide  the  correct 
electrostatic  potential  in  region  I  and  proper  boundary 
conditions  for  the  atoms  at  the  border  of  regions  I  and  II. 
Regions  I  and  II  combined  together  form  a  finite  region 
which  has  a  radius  of  several  nanometers  and  is  called  a 
nano-cluster.  Finally,  to  account  for  the  polarization  of 
the  solid  from  the  border  between  regions  I  and  II  and 
up  to  infinity,  we  introduce  region  III.  It  is  treated  in  the 
approximation  of  a  polarizable  continuum  and  conforms 
geometrically  to  the  boundary  between  regions  I  and  II. 

The  whole  approach  is  similar  in  spirit  to  the  classical 
Mott-Littleton  method47.  The  total  energy  of  the  system 
includes  (see  Appendix  A):  i)  the  energy  of  the  QM  clus¬ 
ter  and  the  interface  atoms  in  the  external  electrostatic 
potentials  due  to  the  classical  environment;  ii)  the  inter¬ 
action  between  classical  atoms  in  regions  I  and  II  calcu¬ 
lated  using  inter-atomic  potentials;  iii)  the  short-range 
interaction  of  the  classical  atoms  and  of  the  interface 
atoms  with  the  QM  cluster  atoms  calculated  using  the 
short-range  part  of  the  classical  inter-atomic  potentials, 
and  iv)  the  Mott-Littleton  correction  for  the  polarization 
of  region  III.  This  latter  correction  is  about  0.3  eV  in  the 
typical  setup  described  below  and  remains  constant  in  all 
calculations  of  singly  charged  vacancies.  The  expression 
for  the  total  energy  is  given  in  ref.31. 

The  classical  atoms  in  regions  I  and  II  are  represented 
using  the  shell  model  and  the  rigid  atom  model,  respec¬ 
tively  (see  Appendix  A).  Their  electrostatic  interaction 
with  the  interface  and  with  the  QM  cluster  atoms  is  in¬ 
cluded  at  the  quantum-mechanical  level,  i.e.  by  calculat¬ 
ing  the  corresponding  matrix  elements  and  adding  them 
to  the  QM  potential  energy  matrix.  Total  forces  are 
calculated  on  both  QM  and  classical  ions.  This  allows 
us  to  minimize  the  system  total  energy  simultaneously 
with  respect  to  the  electronic  coordinates  and  the  posi¬ 
tions  of  QM  ions  and  classical  ions,  and  hence  to  avoid 
the  time-consuming  ’’self-consistency”  procedure  used  in 
some  previous  implementations  of  this  method48,49. 

Properties  of  the  interface  atoms  depend  on  the  ma¬ 
terial  in  question,  as  described  for  MgO39,  Si0231  and 
Mg2SiC>450.  For  Si02,  a  QM  cluster  is  terminated  by  the 
Si*  atoms  that  form  an  interface  between  the  QM  cluster 
and  the  rest  of  the  nano-cluster  (see  Appendix  A).  The 
Si*  atoms  are  located  at  the  Si  sites  of  the  Si02  struc¬ 
ture  and  represent  real  Si  atoms.  This  is  quite  different 
from  many  molecular  cluster  schemes  where  a  QM  clus¬ 
ter  is  terminated  by  artificial  so  called  pseudo-atoms51-53. 
The  Si*  atoms  are  chosen  so  that  they  are  coordinated 
by  one  quantum-mechanically  treated  oxygen  and  three 
classically  treated  oxygen  ions  from  the  rest  of  the  nano¬ 
cluster.  They  perform  a  dual  role:  they  form  a  polar 
bond  with  the  QM  oxygen  and  describe  the  interaction 
with  the  three  classical  oxygens.  This  is  achieved  by 
representing  Si*  as  a  combination  of  j  of  a  quantum- 
mechanical  Si  atom  and  |  of  a  classical  Si  atom.  The 
detailed  description  of  Si*  atoms  as  applied  to  a-quartz 


is  given  in  refs31,32.  The  same  parametrization  scheme  is 
used  in  the  present  calculations. 

The  GUESS  code31,39  plays  the  role  of  a  ’’master”  pro¬ 
gram  that  calculates  the  total  energy,  total  forces  acting 
on  all  centers  in  Region  I  and  performs  the  geometry  op¬ 
timization  of  the  whole  system.  Thus  the  scheme  allows 
us  to  account  consistently  for  the  defect-induced  polar¬ 
ization  of  the  host  lattice  and  also  the  effect  of  this  polar¬ 
ization  on  the  defect  itself.  The  Gaussian  98  package54 
is  used  for  calculations  of  the  quantum  mechanical  con¬ 
tributions  to  the  total  energy  and  forces. 


A.  Summary  of  the  new  method 

The  main  purpose  of  this  work  has  been  to  develop  fur¬ 
ther  the  embedded  cluster  method  for  silica  proposed  in 
refs.31,32’43,  and  to  test  it  on  well-studied  defects.  We 
aim  to  minimize  the  lattice  distortions  at  the  border 
of  a  defect-free  QM  cluster  embedded  in  the  classically 
treated  non-defective  crystal  lattice  and  to  represent  the 
quartz  electronic  structure  as  closely  as  possible  to  that 
calculated  for  the  perfect  lattice  treated  periodically.  For 
future  calculations  of  spectroscopic  properties,  one  would 
like  to  avoid  having  spurious  electronic  states  due  to  the 
cluster  border  at  the  top  of  the  valence  band  and  at  the 
bottom  of  the  conduction  band.  The  scheme  developed  in 
this  project  employs  a  hybrid  density  functional  similar 
to  the  one  known  as  B3LYP  which  includes  the  three- 
parameter  Becke’s  exchange  functional55  and  a  correla¬ 
tion  functional  by  Lee,  Yang,  and  Parr56,  rather  than 
the  Hartree-Fock  method  used  in  the  earlier  studies31,32. 
It  makes  use  of  the  new  semi-local  embedding  potentials 
and  benefits  from  improved  classical  inter-atomic  poten¬ 
tials.  The  latter  provide  a  fully  consistent  treatment 
of  the  crystal  surrounding  the  quantum  cluster  both  in 
terms  of  effective  ionic  charges  and  the  elastic  and  dielec¬ 
tric  properties. 

The  new  generation  of  the  computational  embedding 
technique  developed  within  this  project  includes  several 
components  outlined  below  and  discussed  in  detail  in  Ap¬ 
pendix  A.  They  are  generic  to  further  applications  of  this 
method  to  other  systems,  such  as  chalcogenates  and  high 
dielectric  constant  oxides,  e.g.  Hf02  and  hafnium  sili¬ 
cates.  The  results  of  test  calculations  carried  out  for  the 
non-defective  a-quartz  lattice  and  for  the  structure  and 
spectroscopic  properties  of  the  neutral  oxygen  vacancy  in 
a-quartz  are  summarized  in  Appendix  B.  The  results  of 
applying  of  the  new  method  to  studying  defect  properties 
in  amorphous  silica  are  presented  in  refs.31’32’43’57-60. 

The  computational  scheme  includes  the  following  com¬ 
ponents: 

First,  a  periodic  calculation  of  the  perfect  crystal, 
a-quartz,  is  carried  out  using  the  CRYSTAL  code61. 
This  code  employs  the  Gaussian-type  basis  set  and  im¬ 
plements  the  Hartree-Fock  method  and  hybrid  density 
functionals  compatible  with  those  implemented  in  the 
Gaussian  package  for  quantum-chemical  calculations  of 
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molecules54  and  used  in  further  embedded  cluster  cal¬ 
culations.  In  periodic  calculations  we  used  the  6-31G* 
basis  set  and  optimized  the  geometry  of  the  a-quartz 
lattice  including  both  fractional  coordinates  and  the  unit 
cell  parameters. 

Second,  the  parameters  of  the  classical  inter-atomic 
potentials  are  re-adjusted  to  reproduce  the  structure  ob¬ 
tained  from  the  periodic  calculations  and  experimental 
physical  properties  such  as  elastic  and  dielectric  con¬ 
stants.  The  latter  are  particularly  important  for  accu¬ 
rate  calculations  of  charged  defects.  The  fitted  parame¬ 
ters  provide  consistency  of  the  geometry  in  the  classical 
region  with  that  within  the  QM  cluster.  Changing  basis 
set  or  functional  leads  to  different  geometric  parameters 
of  the  perfect  lattice  and  may  result  in  artificial  strain 
at  the  cluster  border  if  the  classical  interactions  are  not 
adjusted  accordingly.  To  efficiently  parametrize  the  clas¬ 
sical  potentials  for  the  different  basis  set  and/or  ab  initio 
methods  used  in  the  QM  calculations,  we  developed  a 
new  computational  method  which  is  described  in  detail 
in  Appendix  A.  This  approach  allows  us  to  generate  clas¬ 
sical  potentials  that  provide  an  atomic  structure  consis¬ 
tent  with  the  geometry  of  the  QM  cluster  and,  through 
this,  to  minimize  the  lattice  distortions  at  the  interface 
region  in  the  non-defective  system. 

Third,  the  parameters  of  the  interface  atoms,  Si*,  at 
the  border  between  the  QM  cluster  and  the  classical  re¬ 
gion  are  optimized.  These  include  the  parameters  of  a 
semi-local  pseudo-potential  describing  the  interaction  of 
the  border  Si*  atoms  with  the  electron  density  in  the 
QM  cluster  and,  effectively,  saturating  the  Si*-0  bond 
directed  inside  the  cluster.  This  is  the  vital  part  of  the 
embedding  scheme  as  it  is  responsible  for  eliminating 
spurious  border  electronic  states  and  achieving  a  homo¬ 
geneous  electron  density  distribution  inside  the  cluster. 
In  addition,  there  are  short-range  classical  potentials  as¬ 
sociated  with  the  Si*  and  its  oxygen  neighbors.  These 
potentials  aim  to  correct  the  local  atomic  structure  at 
the  interface.  This  approach  allows  us  to  significantly 
improve  the  description  of  the  cluster  border  (compare, 
for  example,  with  previous  publications  by  Sulimov  et 
al.31  and  Mysovsky  et  al.32).  In  particular,  since  the  pa¬ 
rameters  of  the  classical  inter-atomic  potentials  and  the 
effective  pseudo-potential  for  Si*  are  consistent  with  the 
geometry  in  the  QM  region,  the  mismatch  in  the  local 
atomic  structure  at  the  interface  is  minimized. 

Yet  another  difference  between  this  work  and  the  pre¬ 
vious  publications  has  been  that  we  employed  the  DFT 
method  and  a  hybrid  density  functional  similar  to  the 
B3LYP  functional54  instead  of  the  Hartree-Fock  method 
used  before.  We  adjusted  the  amount  of  exact  exchange 
interaction  included  in  this  hybrid  functional  to  consider 
self-trapped  holes  and  to  calculate  optical  absorption  en¬ 
ergies  of  point  defects  in  silica  in  agreement  with  exper¬ 
iment  (see  Appendix  A).  In  combination  with  a  Time- 
dependent  DFT  scheme  for  calculating  optical  transition 
energies,  this  led  to  fairly  accurate  predictions  of  the  op¬ 
tical  transition  energies  for  both  neutral  and  positively 


charged  vacancies  as  demonstrated  in59>60’62.  This  is  a 
considerable  improvement  with  respect  to  the  previous 
studies  which,  paves  the  way  to  more  extensive  studies 
of  optical  properties  of  defects  in  amorphous  silica. 

To  summarize,  the  developed  embedded  cluster 
method  allows  one:  i)  To  study  the  full  extent  of  the 
defect-induced  lattice  distortion;  ii)  To  consider  a  fairly 
large  region  of  the  amorphous  structure  and  thus  to  ex¬ 
plore  a  greater  variety  of  different  sites  in  the  amorphous 
structure  and  to  build  more  complete  statistics  of  differ¬ 
ent  defect  types  and  configurations,  iii)  To  calculate  a 
wide  range  of  defect  properties,  including  the  optical  ab¬ 
sorption  spectra.  The  experience  gained  in  these  calcula¬ 
tions  and  in  combining  quantum  mechanical  and  classical 
techniques  will  be  useful  for  further  studies  of  defects  in 
disordered  materials. 


V.  APPENDIX  A:  THE  EMBEDDED  CLUSTER 
METHOD 

In  this  Appendix  we  describe  in  detail  the  embedding 
method,  the  procedure  for  optimizing  the  classical  shell 
model  inter-atomic  potentials  and  the  embedding  poten¬ 
tial  parameters. 


A.  Calculation  of  the  total  energy 


In  our  method  the  whole  system  is  represented  using 
a  large  finite  nano-cluster  of  either  crystalline  or  amor¬ 
phous  silica.  It  uses  Si(|0)4  structural  elements  as  build¬ 
ing  blocks,  as  described  in  ref.43.  The  nano-cluster  is 
divided  into  regions  I  and  II  as  shown  in  Fig.  2a.  The 
expression  for  the  total  energy  is  given  by: 


E  =  Eqm+if  +  Wqm+if,cl  +  Eql  +  Eml ,  (1) 

where  Eqm+if  is  the  total  energy  of  the  QM  cluster  and 
the  interface  in  the  electrostatic  potential  V§£ul  due  to 
the  classical  environment: 
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Eqm+if  also  includes  a  classical  term  WffqM,  which 
describes  a  non-Coulombic  interaction  between  the  in¬ 
terface  and  QM  cluster  atoms;  it  is  a  correction  to  the 
electronic  pseudo-potentials  Vj  centred  on  the  interface 
atoms.  Similarly,  the  classical  term  Wqm+if  gl  de¬ 
scribes  the  non-Coulombic  short-range  interaction  be¬ 
tween  the  classical  atoms  and  atoms  of  the  QM  cluster 
and  the  interface.  Ecl  is  the  total  energy  of  the  classical 
environment  and  Eml  is  the  Mott-Littleton  polarization 
energy  of  the  lattice  outside  region  I. 


a)  b) 


FIG.  2:  a)  General  setup  for  the  embedded  cluster  calculations:  the  division  of  the  whole  system  into  sub-regions  is  shown,  b) 
Schematic  description  of  the  interface  between  the  quantum-mechanically  treated  cluster  and  classical  environment. 


In  more  detail: 


Eqm+if  =  {mo  +  WS°luqM+IF  l*>  +  W?f,rQM'  (2) 
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The  total  energy  of  the  classical  region  is  given  by 
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In  the  expressions  above  Nqm+if  is  the  total  number 
of  atoms  in  both  the  QM  cluster  and  at  the  interface  to¬ 
gether,  whilst  Ne  stands  for  the  total  number  of  electrons 
shared  between  them.  Ncl  is  the  total  number  of  clas¬ 
sical  atoms,  and  Ncor  and  Nsfle  are  the  total  numbers  of 
classical  cores  and  shells  respectively.  Z  are  the  charges 
of  nuclei  in  the  QM  cluster  and  at  the  interface,  Qcor  and 
Qshe  are  charges  of  the  classical  cores  and  shells.  Rcor 
and  Rs?le  are  the  coordinates  of  the  classical  cores  and 
shells.  The  Hamiltonian  Hq  includes  the  kinetic  energy 
of  the  electrons,  the  electrostatic  interaction  of  the  elec¬ 
trons  and  the  nuclei  in  the  QM  cluster  and  the  interface, 
and  effective  pseudo-potentials  for  the  interface  atoms. 
Terms  WcTJqm+if  and  R if  qm  can  represented  by 
pair-wise  Buckingham  or  Morse-type  classical  potentials 
Uij,  defined  for  each  pair  of  atoms  i  and  j,  as  well  as  3- 
and  4-body  classical  potentials.  Constants  k  are  related 
to  polarizability  of  ions.  The  Mott-Littleton  correction 
Eml  is  calculated  over  a  region  which  conforms  geomet¬ 
rically  to  the  boundary  of  region  I  and  extends  to  infinity. 
The  constants  M represent  on-site  polarizabilities,  Q 
is  the  net  charge  of  the  defect  with  respect  of  the  lattice, 
and  r7;  is  the  distance  from  each  atoms  to  the  centre  of 
region  I63. 


Nqm+if  Ncl 

Wcl,qm+if  =  E  E  ^b(R'  “  Rj)> 

i  j 


The  quantum-mechanical  contributions  to  the  total  en¬ 
ergy  and  forces  are  calculated  using  the  Gaussian  98 
package54  for  ab  initio  calculations  of  molecules.  Details 
of  the  parametrization  of  the  interface  Si*  for  silica  are 
considered  in  more  detail  below  in  section  VB3. 


and 
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B.  Optimizing  the  embedding  parameters 

The  scheme  described  above  requires  an  adequate 
choice  of  the  ab  initio  method  and  basis  set  for  the  QM 
cluster  and  of  the  inter-atomic  potentials  for  the  classi¬ 
cal  atoms  in  regions  I  and  II.  These  should  be  consistent 
with  the  properties  of  the  interface  atoms  at  the  border 
between  the  QM  cluster  and  the  classical  region.  In  this 
section  we  describe  our  approach  to  the  parametrizing 
these  components  of  the  scheme. 


1.  Classical  potentials 

Several  sets  of  classical  inter-atomic  potentials  for  sil¬ 
ica  have  been  previously  developed.  However,  these 
potentials  do  not  fully  satisfy  our  needs  for  two  rea¬ 
sons.  Firstly,  we  wish  to  achieve  consistency  between 
the  atomic  structure  of  the  QM  cluster  and  that  of  the 
classical  region.  In  other  words,  the  classical  potentials 
should  provide  the  same  cell  parameters  and  fractional 
coordinates  as  those  obtained  from  a  reference  quantum- 
mechanical  calculation.  In  our  case  the  reference  cal¬ 
culation  is  carried  out  using  the  CRYSTAL  code  and 
Gaussian-type  basis  set.  Secondly,  to  the  best  of  our 
knowledge,  classical  potentials  with  non-formal  charges 
on  Si  and  O  ions64,65  fail  to  reproduce  the  static  dielec¬ 
tric  constant  of  a-quartz  and  that  of  amorphous  SiC>2. 

Since  the  atomic  structure  obtained  in  these  calcula¬ 
tions  will  depend  on  the  basis  set  and  the  ab  initio  tech¬ 
nique,  a  method  for  fast  re-optimization  of  the  classical 
potentials  is  needed.  This  method  has  been  developed 
on  the  basis  of  an  existing  package  MERLIN66  for  mini¬ 
mization  of  functionals  in  multi-dimensional  space. 

In  the  following  we  describe  the  details  of  the  reference 
CRYSTAL  calculation  used  in  this  work,  the  details  of 
the  re-optimization  procedure,  and,  finally,  compare  the 
performance  of  the  classical  potentials  found  in  this  work 
against  those  reported  previously. 

The  primitive  unit  cell  of  a-quartz  is  defined  by  six 
parameters:  two  of  them  (a  and  c)  determine  the  size 
of  the  unit  cell  and  the  remaining  four  (u,  x,  y  and  z) 
determine  the  fractional  coordinates  of  Si  and  O  atoms. 
The  total  energy  of  a-quartz  was  minimized  with  respect 
to  these  six  parameters  using  the  periodic  model  and  the 
CRYSTAL  code.  Standard  6-31G*  basis  sets  on  Si  and  O 
atoms  and  the  hybrid  B3LYP  density  functional61  were 
employed  in  these  calculations.  The  energy  minimization 
was  done  in  two  steps.  First,  the  fractional  coordinates 
were  fully  relaxed  for  the  fixed  lattice  parameters  a  and 
c.  Then  the  fractional  coordinates  were  fixed  and  a  and 
c  were  relaxed  manually.  Then  both  steps  were  repeated 
until  overall  convergence  was  achieved.  The  total  energy 
per  unit  cell  was  converged  within  ICC6  eV,  the  lattice 
cell  vectors  within  0.0005  A  and  the  fractional  coordi¬ 
nates  within  0.0005.  The  structure  of  a-quartz  obtained 
here  is  in  agreement  with  the  results  of  previous  calcu¬ 
lations67.  In  particular,  we  found  that  a=4.933  A  and 


TABLE  I:  Parameters  of  the  classical  inter-atomic  potential 
optimized  for  a-quartz.  A,  p  and  C  are  the  Buckingham  po¬ 
tential  parameters  for  the  interaction  between  Si-O  and  0-0 
ions,  k  is  the  spring  constant  for  the  interaction  between 
O  core  and  shell.  The  Si  charge  is  2.4|e|,  O  core  charge  is 
—  1.619|e|,  and  O  shell  charge  is  0.419|e|. 


A  (eV) 

P  (A) 

C  (eV-A6) 

k  (eV/A2) 

Si-O 

3185.06 

0.24102 

46.72 

0-0 

7682.67 

0.25627 

0.01 

16.6 

c=5.402  A. 

To  optimize  the  parameters  of  the  inter-atomic  poten¬ 
tials,  p  =  (p\,P2,  pn ),  we  have  defined  a  cost  function 

f°,  f)  =  YlWi'  ~  A(P))2> 

i 

where  f  =  (/i,  /A  ••)  define  numerical  values  of  the  prop¬ 
erties  to  be  fitted  to  the  corresponding  set  f°  of  the  refer¬ 
ence  values.  Weights  Wi  were  set  to  80.0  for  the  geomet¬ 
rical  parameters  of  the  unit  cell  and  to  2.0  for  the  elastic 
and  the  dielectric  constants.  Reference  geometrical  pa¬ 
rameters  in  f°  were  taken  from  the  reference  CRYSTAL 
calculation  while  the  elastic  and  dielectric  constants  were 
taken  from  experimental  data. 

The  optimization  of  parameters  p  was  carried  out  us¬ 
ing  the  MERLIN  optimization  package66  which  imple¬ 
ments  a  modification  of  the  ”  Controlled  Random  Search” 
(CRS)  technique  originally  introduced  by  Price68.  To  cal¬ 
culate  the  value  of  the  cost  function  for  each  probe  set 
of  parameters  p,  we  use  the  GULP  code  (General  Utility 
Lattice  Program)63  which  relaxes  the  a-quartz  lattice  us¬ 
ing  the  inter-atomic  potentials  supplied  by  the  MERLIN 
search  engine  and  then  calculates  the  physical  proper¬ 
ties  at  the  equilibrium  configuration.  After  that  the  cost 
function  for  the  given  probe  set  of  parameters  is  calcu¬ 
lated  and  supplied  back  to  the  search  engine.  In  this 
way,  all  the  parameters  of  the  classical  potentials  for  Si- 
O  and  0-0  interaction,  as  well  as  the  spring  constants 
and  charges  of  the  O  shells  can  be  fitted  simultaneously. 
Thus  this  approach  provides  a  flexible,  yet  efficient  tool 
for  a  general  minimization  problem  able  to  cope  with  very 
complex  landscapes. 

To  fit  the  parameters  for  the  inter-atomic  potentials  we 
first  set  the  total  ionic  charges  on  Si  and  O  ions  to  2.4|e| 
and  — 1.2 |e|  respectively  and  defined  upper  and  lower  lim¬ 
its  for  all  other  parameters.  The  potentials  generated 
using  the  MERLIN  code  are  summarized  in  Table  I. 

The  properties  of  a-quartz  calculated  using  the  pa¬ 
rameters  obtained  in  this  work  are  compared  with  exper¬ 
imental  data  in  Table  II.  The  same  set  of  properties  was 
calculated  using  other  2-body  potentials  reported  previ¬ 
ously.  The  BKS  potentials65  were  modified  so  that  the 
rigid  ion  model  for  oxygen  ions  was  replaced  by  the  po¬ 
larizable  ion  model43.  The  original  potentials  reported 
by  SMC69  also  contain  a  corrective  3-body  term;  the  re- 
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TABLE  II:  Comparison  of  several  parameter  sets  for  a-quartz. 
a  and  c  are  the  hexagonal  lattice  parameters  in  A,  B  is  the 
bulk  modulus  in  GPa  units,  cq  are  the  elastic  tensor  com¬ 
ponents  in  10  GPa  units,  e/,  and  £33  are  the  static  dielectric 
constants,  and  £oo  is  the  high  frequency  dielectric  constant. 


experiment 

this  work 

BKS1 

TTAM2 

SMC3 

a 

4.916 

4.931 

4.933 

5.041 

4.849 

c 

5.405 

5.404 

5.432 

5.552 

5.401 

b 

37 

37.8 

41.8 

32.5 

- 

Cn 

8.69 

8.87 

9.14 

7.25 

6.20 

C13 

1.19 

0.95 

1.73 

1.13 

1.63 

C14 

-1.81 

-1.60 

-1.72 

-1.49 

-1.01 

C33 

10.50 

10.51 

11.13 

8.71 

7.44 

C44 

5.28 

5.63 

5.09 

3.94 

3.30 

Cee 

3.99 

3.90 

4.07 

3.23 

2.74 

e?i 

4.52 

4.41 

1.96 

2.09 

5.51 

e33 

4.46 

4.55 

2.01 

2.15 

6.09 

Coo 

2.4 

2.40 

1.54 

- 

2.06 

Ionic  charges 

Si 

2.4 

2.4 

2.4 

4.0 

O 

-1.2 

-1.2 

-1.2 

-2.0 

1  from  ref.  [65] 

2  from  ref.  [64] 

3  from  ref.  [69] 


suits  presented  in  Table  II  have  been  calculated  without 
the  3-body  correction. 

We  note  that,  as  well  as  a  generally  better  agree¬ 
ment  with  the  experimental  data,  the  potentials  obtained 
in  this  work  reproduce  well  the  reference  geometry,  for 
which  they  were  fitted  and  the  dielectric  constants.  The 
latter,  as  we  discussed  above,  is  essential  for  an  accurate 
calculation  of  the  long-range  lattice  polarization  induced 
by  charged  defects. 


2.  The  density  functional 

It  is  a  well-established  fact  that  different  quantum- 
mechanical  methods  can  produce  qualitatively  different 
results.  For  example,  plane-wave  DFT  calculations  with 
PW91  functional  do  not  reproduce  the  electronic  struc¬ 
ture  of  the  Vk  centre  in  NaCl34,  a  prototype  hole  center 
well-studied  experimentally.  This  has  been  attributed  to 
the  inaccurate  treatment  of  the  exchange  interaction  in 
the  PW91  density  functional70.  On  the  other  hand,  cal¬ 
culations  made  with  the  Hartree-Fock  approach  do  in¬ 
clude  the  exact  exchange  EFXF  but  suffer  from  the  lack 
of  the  electron  correlation.  A  way  forward  was  suggested 
in  the  form  of  hybrid  density  functionals,  which  contain 
an  admixture  of  EFF.  For  example,  there  is  20%  of  EFF 
in  the  B3LYP55  and  50%  in  the  BH&H71  density  func¬ 
tionals.  However,  depending  on  the  system  and/or  class 
of  problems,  the  amount  of  EFF  that  needs  to  be  in¬ 


cluded  in  the  calculations  may  vary.  A  possible  strategy 
is  to  include  the  parametrization  of  the  E^XF  admixture 
that  fits  best  for  a  given  problem. 

We  illustrate  this  idea  with  the  example  of  optical  ab¬ 
sorption  of  the  neutral  vacancy  in  a-quartz.  We  form  a 
BxLYP  density  functional  and  vary  the  E^F  contribu¬ 
tion  so  as  to  reproduce  the  optical  absorption  spectrum 
of  the  S^Hg  molecule72.  This  molecule  has  been  previ¬ 
ously  considered  as  a  model  for  the  neutral  vacancy  in 
a-quartz. 

We  found  that  the  main  features  of  the  S^Hg  absorp¬ 
tion  spectrum  were  reasonably  well  reproduced  if  the 
B3LYP  functional  and  a  standard  6-31G  basis  set  were 
used.  However,  with  a  more  flexible  basis  set  aug-cc- 
pVTZ54,  the  absorption  maxima  shifted  to  lower  ener¬ 
gies.  When  the  amount  of  EFF  was  increased  from  20% 
to  32.5%,  the  agreement  between  the  experimental  and 
theoretical  spectra  improved  if  the  same  aug-cc-pVTZ 
basis  set  was  used.  In  the  following  calculations  (see  Sec¬ 
tion  VI)  we  use  this  modified  BxLYP  functional. 

A  similar  fitting  approach  can  be  used  when,  for  exam¬ 
ple,  properties  of  self-  trapped  holes  in  a-SiC>2  or  excitons 
are  studied.  In  these  cases  the  amount  of  the  EFXF  will, 
generally  speaking,  be  different. 

3.  The  interface  Si*  atoms 

In  the  present  embedded  cluster  calculations  of  SiC>2, 
a  QM  cluster  is  terminated  by  Si*  atoms  that  form  an 
interface  between  the  QM  cluster  and  the  rest  of  the  sys¬ 
tem.  The  Si*  atoms  are  located  at  the  Si  sites  of  the  SiC>2 
structure  and  represent  real  Si  atoms.  This  is  in  contrast 
to  many  other  molecular  cluster  schemes  where  a  QM 
cluster  is  terminated  by  artificial  so-called  pseudo-atoms 
located  in  the  middle  of  Si-0  bonds51-53,73-75. 

The  Si*  atoms  are  chosen  so  that  they  are  coordi¬ 
nated  by  one  quantum-mechanically  treated  oxygen  and 
three  classically  treated  oxygen  ions.  The  dual  role  of 
Si*  atoms  is  to  interact  quantum-mechanically  with  their 
QM  neighbors,  i.e.  to  form  a  polar  bond  with  QM  oxy¬ 
gen  atoms,  and  interact  classically  with  classical  oxygens. 
This  is  achieved  by  representing  Si*  as  a  combination  of 
|  of  a  quantum-mechanical  Si  atom  and  |  of  a  classical 
Si  atom. 

The  detailed  description  of  such  a  division  is  given  in 
refs.  31,32.  Briefly,  Si*  is  an  one-electron  atom;  its  inter¬ 
action  with  the  rest  of  the  system  is  defined  by  several 
parameters  which  are:  i)  the  nucleus  charge  of  Si*  re¬ 
lated  to  the  effective  ionic  charge  of  silicon  in  SiC>2,  ii) 
the  effective  core  pseudo-potential  V)  (r),  iii)  the  atomic 
basis  set,  and  iv)  the  classical  short-range  inter-atomic 
potentials  that  describe  the  interaction  of  the  Si*  with  its 
neighbors  in  the  QM  cluster  and  in  the  classical  environ¬ 
ment.  These  components  are  inter-related  and,  therefore, 
should  be  fitted  together. 

The  meaning  of  the  first  two  components  of  the  Si* 
parametrization  can  be  highlighted  as  follows.  For  a 
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given  ionic  charge  Qsi  for  Si  ions  in  SiC>2  the  amount  of 
the  electron  charge  donated  to  each  neighboring  O  is  cal¬ 
culated  as  \Qsi-  Therefore,  Si*,  as  an  one-electron  atom 
with  only  one  quantum-mechanically  treated  neighbor, 
should  have  a  nuclear  charge  equal  to  1  +  '\Qs-i-  This 
provides  a  deep  potential  well  for  the  electrons,  which 
has  to  be  locally  compensated  by  a  short-range  Effective 
Core  Pseudo-potential  (ECP)  to  reproduce  the  effect  of 
the  core  electrons  in  real  Si  atoms.  The  atomic  Gaussian- 
type  basis  set  of  Si*  consists  of  a  single  s  basis  function 
that  corresponds  to  the  valence  shell  of  the  6-31G  basis 
set  for  Si  atoms. 

In  the  current  parametrization  we  employed  a  semi¬ 
local  pseudo-potential,  which  is  better  suited  for  de¬ 
scribing  the  directed  Si-0  bond  than  the  local  pseudo¬ 
potential  used  in  the  previous  parametrization.  The  an¬ 
alytical  form  of  this  pseudo-potential  is  that  of  Hay  and 
Wadt19: 

L—  1 

P(r)  =  UL(r)  +  Y^mr)  -  UL{r)\Pi){F\\, 

1=0 

where  \Pi){Pi\  is  the  projector  on  the  basis  functions  of 
the  orbital  momentum  l. 

In  this  project  we  developed  the  Si*  ECP  with  L  =  1. 
Functions  Ui(r)  were  expanded  over  10  primitive  func¬ 
tions  of  the  form:  Ce~ar  r2~n  where  C  are  expansion 
coefficients  and  n  =  0, 1, 2. 

To  fit  the  potential  we  fixed  the  powers  of  rij  and  var¬ 
ied  parameters  Ct  and  a,  in  two  steps.  First,  the  effec¬ 
tive  potential  for  Si*  was  fitted  on  the  Hartree-Fock  level 
for  the  SQOySig  QM  cluster  in  the  ideal  a-quartz  con¬ 
figuration.  In  these  calculations  we  demanded  that  the 
distribution  of  the  electron  density,  as  given  by  the  Nat¬ 
ural  Population  Analysis  (NPA),  inside  the  QM  cluster 
was  homogeneous  and  the  direction  and  the  electronic 
structure  of  Si*-0  bonds,  as  given  by  Natural  Bonding 
Analysis  (NBA),  were  similar  to  those  calculated  for  real 
Si-0  bonds.  Moreover,  the  average  charges  of  the  three 
types  of  ions  should  be  balanced.  This  can  be  expressed 
by  the  following  equations: 

Qsi  =  -2  Qo  =  4  Qsi* 

In  addition,  we  demand  that  the  energy  of  the  localized 
electronic  hole  ( h+ )  should  be  independent  of  whether 
the  hole  is  localized  on  the  oxygen  ion  at  the  centre  of 
the  QM  cluster  or  at  its  boundary.  The  energy  of  the 
hole  was  calculated  as  the  difference  between  the  total 
energies  of  the  ground  and  ionized  states.  A  similar  cri¬ 
terion  was  used  for  the  localized  unrelaxed  excitons  where 
we  demanded  that  the  excitation  energies,  calculated  via 
Singlet-Triplet  transitions  (S'— >T),  were  independent  of 
the  exciton  localization  site. 

The  optimized  set  of  parameters  C),  a*  and  n,  provided 
E{htentre)  =14.93  eV  and  E{h+eriphery)  =14.92  eV.  Simi¬ 
larly,  E(S  ^T)cen£re  =  10.06  eV  and  E(S  *T) periphery  — 


9.57  eV  were  calculated.  These  results  indicate  that  the 
effective  pseudo-potential  developed  for  Si*  in  this  work 
provides  a  reasonably  accurate  embedding  potential  with 
an  expected  error  of  about  0.5  eV  for  the  excitation  en¬ 
ergies  close  to  the  inter-band  transitions. 

In  the  second  step  of  the  Si*  ECP  optimization,  the 
radial  functions  U[  (r)  obtained  above  were  scaled  with 
coefficients  (3  and  7  so  that  C[  =  /3Ct  and  a'  =  7a j, 
where  (3  and  7  were  varied  independently  for  each  Ui(r) 
to  satisfy  the  same  criterion  on  the  homogeneous  distri¬ 
bution  of  the  electron  density  as  described  above.  The 
geometrical  parameters,  e.g.  Si-0  inter-atomic  distances, 
inside  the  QM  cluster  should  were  also  required  not  to 
depend  on  the  distance  from  the  center  of  the  bond  to 
the  boundary  of  the  cluster.  To  achieve  this,  parameters 
f3  and  7  for  each  Ui  (r)  were  scaled  independently  to  mini¬ 
mize  cost-functions  associated  with  the  above  conditions. 
The  MERLIN  package  was  used  for  these  calculations  in 
the  same  spirit  as  described  above  in  section  VB  1. 

Finally,  the  classical  interaction  of  the  interface  Si* 
atoms  with  other  lattice  atoms  was  calculated  as  follows. 
The  pair-wise  interaction  of  Si*  with  classical  atoms  was 
taken  to  be  the  same  as  the  Si-0  interaction  in  the  clas¬ 
sically  treated  part  of  the  system.  The  Si*-0  bond  with 
the  QM  oxygen  was  simulated  using  a  Morse-type  poten¬ 
tial  V(r)  =  De  [(l  —  e-Q(r-'ro))“  —  ll  with  parameters 


De= 28  eV,  a=3.9  A-1,  and  ro=1.59  A.  In  addition  to 
the  pair-wise  interactions,  classical  three-body  potentials 
of  the  form  |(0— 90)2  were  introduced  at  the  interface  be¬ 
tween  the  QM  cluster  and  the  classical  region  to  correct 
for  the  distortion  of  Si*-0-Si  angles.  The  force  constant 
of  k= 4  eV/rad2  was  fitted  for  these  potentials. 


At  the  end  of  the  total  energy  minimization  with  re¬ 
spect  to  the  coordinates  of  all  centers  in  region  I,  we 
found  that  the  local  geometry  of  the  QM  cluster  was  very 
close  to  the  ideal  a-quartz  structure  (see  Table  III) .  The 
largest  deviation  of  an  O  atomic  charge  from  the  corre¬ 
sponding  average  was  smaller  than  0.02  |e|  and  the  ratio 
of  calculated  average  atomic  charges  on  Si,  O,  and  Si* 
atoms  was  3.98:-1.98:l,  i.e.  very  close  to  the  ideal  ratio 
of  4:-2:l.  Details  of  the  perfect  lattice  test  calculations 
for  QM  clusters  SQOySig  and  Si8025Sif8  are  described 
in  the  next  section. 


VI.  APPENDIX  B:  TEST  CALCULATIONS 

To  test  the  performance  of  the  new  method  and  to  il¬ 
lustrate  the  significance  of  different  components  of  the 
embedding  potential,  we  consider  two  types  of  problems. 
First,  to  demonstrate  improved  consistency  between  the 
quantum-mechanically  and  classically  treated  regions,  we 
study  the  electronic  and  geometric  structures  of  the  per¬ 
fect  a-quartz  lattice  as  a  function  of  the  cluster  size  and 
compare  the  results  with  the  periodic  CRYSTAL  cal¬ 
culations  as  well  as  with  the  results  obtained  using  an 
earlier  parametrization  for  Si02  based  on  a  local  effec- 
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TABLE  III:  Average  geometrical  parameters  in  the  QM  cluster  and  at  the  interface  (Si-0  and  Si*-0  distances,  Si-O-Si  and 
Si*-0-Si  angles)  and  NPA  atomic  charges  used  to  characterise  the  quality  of  embedding  in  I4_zoc-scheme,  TA-ioc-scheme, 
and  H-scheme.  Ratio  Q(Si):Q(0,):Q(Si*)  is  calculated  with  respect  to  the  Q(Si*).  Egap  is  the  difference  between  the  highest 
occupied  and  the  lowest  unoccupied  one-electron  states.  The  notations  are  described  in  the  text.  Distances,  angles,  and  the 
atomic  charges  are  given  in  A,  degrees  and  |e|  respectively.  Egap  is  given  in  eV.  Numbers  in  brackets  indicate  the  range  of  the 
corresponding  values  (from  their  maximum  to  minimum). 


Property 

This  work 

Earlier  scheme31 

Free  clusters 

Si207Sig 

SigOasSilg 

Si207Si£ 

SigOasSilg 

Si207H6 

Si8026Hi8 

Si-O 

1.635 

(0.022) 

1.638 

(0.042) 

1.608 

(0.043) 

1.617 

(0.068) 

1.688 

(0.015) 

1.680 

(0.034) 

o 

1 

m 

1.636 

(0.002) 

1.636 

(0.002) 

1.554 

(0.002) 

1.554 

(0.002) 

0.966 

(0.008) 

0.965 

(0.013) 

Si-O-Si 

146.1 

(n/a) 

145.0 

(3.9) 

147.7 

(n/a) 

148.1 

(6.2) 

146.8 

(n/a) 

148.3 

(12.2) 

Si*-0-Si 

142.5 

(1.9) 

142.3 

(2.6) 

149.1 

(1.6) 

147.8 

(3.6) 

131.1 

(11.6) 

128.7 

(16.4) 

Q(Si) 

2.51 

(0.00) 

2.52 

(0.03) 

2.61 

(0.01) 

2.59 

(0.04) 

2.42 

(0.01) 

2.45 

(0.12) 

Q(O') 

-1.27 

(n/a) 

-1.27 

(0.01) 

-1.22 

(n/a) 

-1.25 

(0.03) 

-1.28 

(n/a) 

-1.28 

(0.01) 

Q(o") 

-1.25 

(0.01) 

-1.25 

(0.01) 

-1.41 

(0.01) 

-1.41 

(0.01) 

-1.13 

(0.02) 

1.12 

(0.03) 

Q(Si*)  /  Q(H) 

0.63 

(0.01) 

0.63 

(0.01) 

0.74 

(0.01) 

0.74 

(0.01) 

0.53 

(0.01) 

0.53 

(0.02) 

Q(Si:0':Si*) 

3.98: 

-2.02:1 

4.00:- 

-2.02:1 

3.53: 

1.65:1 

3.50:- 

1.69:1 

4.57:- 

-2.42:1 

4.62: 

-2.40:1 

Egap 

9.34 

8.71 

7.50 

7.04 

10.17 

9.44 

five  pseudo-potential31  (V)oc-scheme),  and  a  traditional 
scheme  where  hydrogen  atoms  are  used  to  terminate  the 
dangling  bonds  (H-scheme).  Second,  the  ground  state 
structure  and  optical  properties  of  the  neutral  oxygen 
vacancy  in  a-quartz  are  calculated  using  the  scheme  pre¬ 
sented  above  and  compared  with  the  results  obtained 
both  Vzoc-scheme  and  H-scheme. 


The  general  setup  for  the  embedded  cluster  calcula¬ 
tions  has  been  described  previously  in31,32,43.  Here  we 
use  a  spherical  nano-cluster  of  radius  of  30  A  built  of 
Si(|0)4  building  blocks.  Region  I,  where  all  atoms  were 
fully  relaxed,  contained  the  same  number  of  atoms  (789) 
in  both  calculations  of  the  non-defective  lattice  and  of 
oxygen  vacancies.  In  the  free  cluster  calculations  the 
terminating  hydrogen  atoms  were  fixed  at  a  distance 
of  0.98  A  from  their  oxygen  neighbours  so  that  the  di¬ 
rections  of  O-H  bonds  coincided  with  the  directions  of 
the  equivalent  O-Si  bonds  in  a-quartz  and  all  other 
atoms  were  fully  relaxed.  Three  QM  clusters  were  used: 
Si207Sig,  Sig025Si*8,  and  S^iOseSi^g,  where  Si*  is  an  in¬ 
terface  atom  in  the  case  of  embedded  cluster  calculations 
or  a  hydrogen  atom  in  the  case  of  the  H-scheme. 

The  standard  6-3 1G  basis  set  was  used  for  all  QM 
atoms  unless  otherwise  stated.  The  hybrid  BxLYP 
functional  containing  the  Lee,  Yang  and  Parr  corre¬ 
lation56  and  modified  Becke  three-parameter  exchange 
functional55  with  32.5%  Hartree-Fock  exchange  was  used. 
The  optical  absorption  energies  were  calculated  using  the 
Time-Dependent  DFT  (TDDFT)  as  implemented  in  the 
Gaussian  98  code54 


A.  Perfect  lattice 

To  assess  the  quality  of  the  embedded  cluster  scheme 
we  have  measured  several  characteristics  listed  in  Ta¬ 
ble  III.  The  geometry  of  the  QM  cluster  after  the  full 
energy  minimization  can  be  described  using:  i)  average 
Si-0  bond  length  and  average  Si-O-Si  angle  inside  the 
QM  cluster,  ii)  average  Si*-0  bond  length  at  the  inter¬ 
face,  and  iii)  average  Si*-0-Si  angle  at  the  QM  cluster 
border.  To  characterise  the  distribution  of  the  electron 
density  inside  the  QM  cluster  one  can  use  NPA  atomic 
charges  averaged  over:  i)  all  Si  atoms  in  the  QM  cluster 
(Q(Si)),  ii)  O  atoms  (denoted  as  O')  that  are  coordinated 
by  two  QM  Si  atoms  in  the  QM  cluster  (Q(O')),  iii)  O 
atoms  (denoted  as  O")  that  are  coordinated  by  one  Si* 
and  QM  Si  atoms  (Q(0")),  and  iv)  all  Si*  atoms  (Q(Si*)). 
Finally,  one  can  compare  the  ratio  of  the  atomic  charges 
Q(Si):Q(0'):Q(Si*)  with  the  ’’ideal”  ratio  of  4  :  — 2  :  1. 

A  comparative  analysis  of  the  data  collected  in  Ta¬ 
ble  III  suggests  that  the  scheme  presented  in  this  work 
provides  a  substantial  improvement  over  the  earlier 
scheme  reported  in  ref.  31.  In  particular,  the  mismatch 
between  Si-0  and  Si*-0  inter-atomic  distances  has  been 
removed  and  the  spread  of  the  geometrical  parameters 
has  become  smaller.  More  importantly,  however,  the 
distribution  of  the  electron  density,  as  manifested  by 
the  atomic  charges,  is  more  homogeneous  in  the  present 
scheme  and  the  ratio  of  charges  is  closer  to  the  stoi¬ 
chiometry  of  the  system.  For  example,  atomic  charges 
on  oxygens  coordinated  by  two  QM  Si  atoms  and  by  a 
QM  Si  and  interface  Si*  atoms  are  almost  identical,  in 
contrast  to  the  earlier  V)oc  scheme.  In  addition,  the  one- 
electron  band  gap  is  now  closer  to  the  experimental  value 
of  9.65  eV  as  determined  by  EELS'7.  The  traditional  H- 
scheme  results  in  a  much  less  accurate  atomic  structure 
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FIG.  3:  Densities  of  states  calculated  for  a-quartz  structure  using  the  embedded  cluster  and  periodic  approaches,  a)  embedded 
cluster  DOS  calculated  using  6-31G  basis  set76  (basG9 s);  b)  periodic  model  DOS  calculated  using  the  same  basis  basG98 ;  c) 
periodic  model  DOS  calculated  using  a  more  localized  basis  developed  for  CRYSTAL  calculations61  ( basco3 ) 


with  the  spread  of  Si-O-Si  angles  as  large  as  12°  and  an 
inhomogeneous  charge  density. 

To  demonstrate  the  dependence  of  the  density  of  states 
(DOS)  on  the  cluster  size  we  compare  the  DOS  calcu¬ 
lated  using  the  current  V's_;oc-scheme  for  two  QM  clusters 
(Fig.  3a)  and  that  calculated  using  the  periodic  model 
and  the  CRYSTAL-03  code  (Figs.  3b  and  3c).  For  the 
sake  of  comparison  we  used  the  same  geometrical  param¬ 
eters  in  all  three  calculations. 

The  DOS  in  Fig.  3a  was  calculated  using  the  embed¬ 
ded  cluster  approach  for  QM  clusters  Sis025Si);8  (black) 
and  Si2i056Si28  (red)  and  a  standard  6-31G  basis  set  for 
molecular  calculations  ( basam  thereafter)76.  The  DOS 
in  Fig.  3b  was  calculated  using  the  periodic  model  and 
the  CRYSTAL-03  code  and  the  same  basc9 s  basis  set. 
Finally,  Fig.  3c  shows  the  DOS  calculated  with  the  same 
periodic  model  and  the  basis  optimized  for  bulk  calcu¬ 
lations  of  Si02  crystal  (base 03)-  Vertical  dashed  lines 
indicate  the  boundaries  of  the  valence  and  conduction 
bands.  The  energy  axes  were  shifted  so  that  the  top  of 
the  valence  band  coincides  in  all  three  cases. 

The  positions  of  the  top  of  the  valence  band  and  the 
bandgap  obtained  in  the  two  embedded  cluster  calcula¬ 
tions  are  identical  in  both  calculations.  The  profiles  of 
the  upper  part  of  the  valence  band  (between  -17  eV  and 
-13  eV  in  Fig.  3a  are  also  very  similar.  In  fact,  the  only 
difference  appears  to  be  in  the  profiles  of  the  lower  part 
of  the  valence  band  (from  -  26  eV  to  -20  eV)  formed  by 


the  Si-0  bonding  states.  Analysis  of  these  states  sug¬ 
gests  that  the  state  associated  with  Si*-0  bonds  at  the 
interface  are  shifted  to  lower  energies  (from  -26  to  -23 
eV)  than  the  equivalent  states  at  the  center  of  the  cluster 
(from  -23  to  -19  eV).  This  shift  is  due  to  the  semi-local 
effective  potential  defined  for  interface  Si*  atoms  which 
’’drives”  the  VB  states  to  lower  energies.  Since  the  rel¬ 
ative  number  of  interface  Si*-0  in  S^iOsgSi^g  cluster  is 
smaller  than  in  Sig025SiJ8  cluster,  the  DOS  in  the  lower 
energy  part  of  the  valence  band  for  the  former  is  more 
accurate  and  shows  better  agreement  with  the  DOS  ob¬ 
tained  in  periodic  calculations.  Finally  we  note  the  re¬ 
sults  of  the  periodic  calculations  are  virtually  identical 
for  the  basis  sets  base 98  and  basco3- 


We  conclude  that  the  band  gap  and  the  main  features 
of  the  Density  of  State  obtained  in  the  embedded  cluster 
calculations  agree  well  with  the  DOS  obtained  in  peri¬ 
odic  calculation  with  the  exception  of  the  Si*-0  bond 
states  at  the  cluster  interfaces.  The  profile  of  the  upper 
VB  remains  almost  unaffected  by  the  cluster  size  thus 
indicating  that  the  nature  of  the  lone  pairs  is  well  repro¬ 
duced  in  relatively  small  clusters.  In  other  words,  the 
main  features  of  the  density  of  states  are  determined  by 
the  interactions  between  the  nearest  Si  and  O  ions  how¬ 
ever  the  valence  band  width  and  fine  features  of  the  DOS 
depend  on  the  cluster  size  and  the  quality  of  embedding. 
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TABLE  IV:  Properties  of  the  neutral  oxygen  vacancy  in  a-quartz.  Average  Si-Si  distances,  position  of  the  vacancy  level  above 
the  Valence  Band  for  the  unrelaxed  (u)  and  relaxed  (r)  vacancies,  unrelaxed  and  relaxed  vacancy  formation  energies  and  optical 
absorption  energies  and  oscillator  strengths. 


Properties 

This  work 

Earlier  scheme31 

Free  clusters 

Si207Sig 

SisCHsSijg 

SbOrSig 

SisCHsSijg 

Si207H6 

SisOasHig 

Ground  state 

Si-Si 

(r) 

2.375 

2.404 

2.383 

2.398 

2.568 

2.398 

evac  —  f-VB 

(u) 

2.55 

2.47 

2.60 

2.84 

2.67 

2.10 

^■vac  B 

(r) 

0.60 

0.49 

0.68 

0.63 

0.72 

0.34 

€Vac  Shift 

1.95 

1.98 

1.92 

2.21 

1.95 

1.76 

E  form 

(u) 

7.96 

7.81 

7.94 

8.04 

7.82 

7.83 

E  form 

(r) 

6.20 

6.18 

6.33 

6.27 

6.78 

6.34 

Erelax 

1.76 

1.63 

1.61 

1.77 

1.04 

1.49 

Optical  transitions 

a  — ►  a* 

(r) 

7.42(0.19) 

7.59(0.44) 

5.24(0.15) 

6.20(0.25) 

7.15(0.40) 

<7  — >  7T 

M 

7.87(0.08) 

7.70(0.19) 

5.86(0.04) 

7.29(0.08) 

7.36(0.05) 

(7  — ►  7T 

8.05(0.09) 

7.81(0.08) 

6.18(0.06) 

7.73(0.04) 

7.62(0.08) 

B.  Properties  of  the  neutral  oxygen  vacancy  in 
a-quartz 

The  properties  of  the  neutral  oxygen  vacancy  calcu¬ 
lated  using  the  three  computational  schemes  are  summa¬ 
rized  in  Table  IV. 

The  neutral  oxygen  vacancy  (NOV)  in  a-quartz  is 
present  in  the  form  of  a  Si-Si  bond  formed  by  two  Si 
atoms  after  an  oxygen  atom,  occupying  the  site  between 
them,  has  been  removed.  It  is  convenient  to  charac¬ 
terize  the  geometric  structure  of  NOV  using  the  dis¬ 
tance  between  the  two  Si  atoms.  Formation  of  the  Si-Si 
bond  induces  a  strong  relaxation  of  the  lattice  around 
the  vacancy31.  Calculations  made  using  both  embedding 
schemes  and  large  H-terminated  clusters  result  in  almost 
identical  Si-Si  bond  lengths.  Clearly,  the  Si-Si  bond  is 
under-relaxed  if  only  a  few  atoms  near  the  vacancy  are 
allowed  to  relax  as  is  the  case  for  small  H-terminated 
clusters. 


An  accurate  account  of  the  NOV-induced  lattice  relax¬ 
ation  is  essential  since  it  affects  the  formation  energy  of  a 
defect  differently  to  its  electronic  properties.  For  exam¬ 
ple,  all  three  schemes  suggest  that  the  un-relaxed  NOV 
formation  energy  is  between  7.8  and  8.0  eV.  After  relax¬ 
ation  both  embedding  schemes  give  the  formation  energy 
as  between  6.2  and  6.3  eV  regardless  of  the  size  of  the  QM 
cluster,  while  in  the  H-terminated  scheme  the  formation 
energy  is  6.3  eV  for  the  larger  cluster  and  6.8  eV  for  the 
smaller  one.  The  calculated  relaxation  energies  suggest 
that  even  in  the  case  of  the  large  H-terminated  QM  clus¬ 
ter  the  NOV  is  not  fully  relaxed,  since  the  corresponding 
relaxation  energy  of  1.5  eV  is  slightly  larger  then  the 
1.6-1. 8  eV  relaxation  for  the  embedding  schemes.  At  the 
same  time,  all  the  computation  schemes  considered  agree 


on  the  position  of  the  one-electron  level  evac  of  the  un¬ 
relaxed  vacancy  and  put  it  at  about  2. 5-2. 8  eV  above 
the  top  of  the  Valence  Band  eyB-  This  level  shifts  down 
by  aproximately  the  same  value  (1.95-2.2  eV)  once  the 
lattice  relaxation  is  taken  into  account. 

The  optical  transitions  calculated  for  NOV  have  three 
major  contributions.  The  most  intensive  transition  cor¬ 
responds  to  an  excitation  from  the  occupied  a  state  of  the 
vacancy  to  the  unoccupied  a*  state  (see  Fig.  4a).  Two 
other  transitions  can  be  characterized  as  a  — ■>  it  transi¬ 
tions  (Fig.  4b).  The  a*  and  two  tt  states  are  induced  by 
the  vacancy  and  their  one-electron  energies  split  slightly 
from  the  bottom  of  the  Conduction  Band.  The  full  op¬ 
tical  absorption  spectrum  also  contains  transitions  from 
the  Valence  Band  states  to  the  unoccupied  states  induced 
by  the  vacancy  and  from  the  occupied  cr  state  to  the  Con¬ 
duction  Band.  These  transitions,  however,  have  much 
smaller  oscillator  strengths. 

The  excitation  energies  and  the  oscillator  strengths  cal¬ 
culated  for  the  most  intensive  transitions  are  summarized 
at  the  end  of  Table  IV.  It  is  immediately  clear  that 
the  excitation  energies  are  greatly  underestimated  if  the 
calculations  are  carried  out  using  the  Si*  parametriza- 
tion  reported  in31.  This  is  not  surprising  since  that 
parametrization  was  based  on  ground  state  calculations 
only.  The  excitation  energies  obtained  using  the  current 
parametrization  are  in  line  with  experimentally  observed 
optical  absorption  of  NOV  at  7.6  eV.  An  extensive  dis¬ 
cussion  on  the  optical  absorption  of  NOVs  in  amorphous 
silica  can  be  found  in60.  Finally,  in  the  case  of  the  small 
H-terminated  free  cluster  S^C^Sig  the  transition  ener¬ 
gies  are  underestimated  with  respect  to  the  experiment. 
As  the  cluster  is  increased,  the  excitation  energies  be¬ 
come  larger  and  closer  to  the  expermental  data.  This 
effect  is  attributed  to  an  insufficiently  accurate  electro¬ 
static  potential  at  the  vacancy  site  and  its  environment 
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FIG.  4:  The  most  active  optical  transitions  for  the  neutral  oxygen  vacancy  in  a-quartz.  a)  cr  — >  a*  transition;  b)  two  u  — >  7r 
transitions. 


in  the  case  of  the  free  clusters.  Clearly,  this  potential  is  taken  into  account  in  the  embedded  cluster  calculations. 
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